---
title: "Replication Data for *Legislative Politics under 'One Country, Two Systems': Evidence from Macau, 2013–2021*"
author: "Daina Chiba"
format:
  html: 
    toc: true
    page-layout: full
    embed-resources: true
editor: visual
date: today
code-fold: false
execute:
  message: false
  warning: false
---

# Preparation

## Packages

Install {pacman} if not installed.

```{r}
if (!requireNamespace("pacman", quietly = TRUE)) {
  install.packages("pacman")
}
```

Load necessary packages.[^1]

[^1]: This replication archive was originally created using ggplot2 version 3.5.2. When running the code with ggplot2 version 4.0.0 or later, the histogram output differs from the original published figure due to changes in the package's binning algorithm or boundary handling between versions. I thus modified the code to produce Figure 7 on 18 October 2025 to ensure that the output remains consistent with the original figure. 

```{r, message = FALSE}
pacman::p_unload(pacman::p_loaded(), character.only = TRUE)
pacman::p_load(tidyverse, showtext, ggpattern, knitr, corrplot, 
               igraph, ggraph, brms, 
               survival, stargazer)

## Display Chinese characters
showtext_auto()
```

## Data

Load legislator-level, bill-level, and vote-level data sets:

```{r, message = FALSE}
df_legislators <- read_csv("Data/df_legislators.csv")
df_bills <- read_csv("Data/df_bills.csv")
df_votes <- read_csv("Data/df_votes.csv")
```

# Figure 1: Legislators

Calculate the number of legislators at each juncture.

```{r}
lgt_totals <- df_legislators |> 
  filter(term != 6 | Legislator_E != "Ho Iat Seng") |> 
  group_by(term, election_type) |> 
  summarize(total = n()) |> 
  filter(term != 4 | election_type != "Direct") |> 
  filter(term != 5 | election_type != "Direct") |> 
  filter(term != 7 | election_type != "Direct") |> 
  filter((term %in% c(2,4))==FALSE | election_type != "Indirect") |> 
  filter((term %in% c(5,7))==FALSE | election_type != "Indirect") |> 
  filter(term == 4 | election_type != "Appointed")
```

Change the ordering of categories:

```{r}
df_legislators <- df_legislators |> 
  mutate(election_type = ordered(election_type, levels = c("Direct", "Indirect", "Appointed")),
         election_type1 = ordered(election_type1, levels = c("D", "I", "A")),
         party = ordered(party, 
                         c("Peak", "Hometown", "Gaming", "Pro-democracy", "Other"))
  )

lgt_totals <- lgt_totals |> 
  mutate(election_type = ordered(election_type, levels = c("Direct", "Indirect", "Appointed"))
         )
```

## Plot

Draw a bar chart that shows the distribution of legislators’ selection methods and political affiliations.

```{r}
#| label: figure-1
#| fig-cap: "Figure 1: Distribution of legislators’ selection methods and political affiliations. This figure shows increases in directly and indirectly elected legislators from 1999 to 2021. The x-axis shows legislative terms since 1999, while the y-axis shows the number of legislators. The distribution changed from 8-8-7 in 1999 to 14-12-7 in 2013."
#| cap-location: margin
#| fig-width: 10
#| fig-height: 6

df_legislators |> 
  filter(term != 6 | Legislator_E != "Ho Iat Seng") |> 
  mutate(term_char = if_else(term == 1, "Term 1: 1999-", NA),
         term_char = if_else(term == 2, "Term 2: 2001-", term_char),
         term_char = if_else(term == 3, "Term 3: 2005-", term_char),
         term_char = if_else(term == 4, "Term 4: 2009-", term_char),
         term_char = if_else(term == 5, "Term 5: 2013-", term_char),
         term_char = if_else(term == 6, "Term 6: 2017-", term_char),
         term_char = if_else(term == 7, "Term 7: 2021-", term_char)
         ) |> 
  ggplot(mapping = aes(x = term_char, fill = party)) + 
  geom_bar_pattern(
    aes(pattern = party),
    pattern_color = "white",
    pattern_fill = "white",
    pattern_spacing = .05, 
    pattern_density = .01
  ) + 
  geom_text(data = lgt_totals, 
            aes(x = term, y = total + .5, label = total, fill = NULL), size = 6) + 
  labs(x = NULL, 
       y = "Number of legislators", 
       fill = "Association", pattern = "Association") + 
  scale_fill_viridis_d() + 
  scale_pattern_manual(values = c("none", "stripe", "circle", "none", "none")) + 
  facet_wrap(election_type ~ .) + 
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 30, vjust = 1, hjust=1)) + 
  theme(legend.key.size = unit(0.8, 'cm'), 
        legend.position = "bottom", 
        strip.text = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        axis.text.x = element_text(size = 10),
        axis.title = element_text(size = 17),
        legend.title = element_text(size = 17),
        legend.text = element_text(size = 15)
        )

ggsave("Figures/Figure1_legislators.pdf", width = 10, height = 6)
```

# Table 1: Roll-Call Record

We will first split the data into General Voting and Detailed Voting parts. We then match the two so that we know the proposers for the Detailed Voting data.

```{r}
df_bills_gen <- df_bills |> 
  filter(vote_type == "General") |> 
  select(BillID, BillTitle, proposer)
df_bills_det <- df_bills |> 
  filter(vote_type == "Detailed") |> 
  select(BillID, BillTitle) |> 
  mutate(x = 1)
temp <- full_join(df_bills_det, df_bills_gen)
temp <- temp |> 
  filter(!is.na(x)) |> 
  select(-x) |> 
  unique() |> 
  select(BillID, BillTitle, proposer_d = proposer)

df_bills <- full_join(df_bills, temp)
df_bills <- df_bills |> 
  mutate(proposer = if_else(is.na(proposer), proposer_d, proposer), 
         vote_resU = if_else(vote_unanimous == TRUE, "Unanimous", "Divided"),
         vote_resP = if_else(vote_pass == TRUE, "Pass", "Fail"), 
         vote_resP = ordered(vote_resP, levels = c("Pass", "Fail")),
         Term_lab = if_else(Term == 5, "Term 5", "Term 6"), 
         Term_lab = ordered(Term_lab, levels = c("Term 6", "Term 5")))
```

## Term 5, General

```{r}
df_bills |> 
  group_by(Term, vote_type, proposer, vote_resU, vote_resP) |> 
  tally() |> 
  filter(Term == 5 & vote_type == "General") |> 
  kable()
```

## Term 6, General

```{r}
df_bills |> 
  group_by(Term, vote_type, proposer, vote_resU, vote_resP) |> 
  tally() |> 
  filter(Term == 6 & vote_type == "General") |> 
  kable()
```

## Term 5, Detailed

```{r}
df_bills |> 
  group_by(Term, vote_type, proposer, vote_resU, vote_resP) |> 
  tally() |> 
  filter(Term == 5 & vote_type == "Detailed") |> 
  kable()
```

## Term 6, Detailed

```{r}
df_bills |> 
  group_by(Term, vote_type, proposer, vote_resU, vote_resP) |> 
  tally() |> 
  filter(Term == 6 & vote_type == "Detailed") |> 
  kable()
```

## Total

```{r}
df_bills |>
  group_by(proposer, vote_type, vote_resU, vote_resP) |>
  tally() |>
  filter(vote_type == "Detailed") |> 
  kable()
```

```{r}
with(df_bills |> 
       filter(vote_type == "General" & Term == 5), table(proposer, vote_resU))
```

```{r}
df_bills |> 
  filter(Term == 6) |>
  filter(proposer == "Government") |>
  filter(vote_type == "Detailed") |>
  pull(vote_resU) |>
  table() |> 
  kable()
```

```{r}
df_bills |> 
  filter(Term == 6) |>
  filter(proposer == "Government") |>
  filter(vote_type == "Detailed") |>
  filter(vote_resU == "Divided") |>
  pull(vote_pass) |>
  table() |> 
  kable()
```

# Figure 2: Legislator-Proposed Bills

Focus on legislator-proposed bills:

```{r}
df_bills_legis <- df_bills |> 
  filter(proposer == "Legislator" & vote_type == "General") |> 
  mutate(proposer_two = if_else(!is.na(Proposer2), TRUE, FALSE),
         proposer_three = if_else(!is.na(Proposer3), TRUE, FALSE),
         proposer_four = if_else(!is.na(Proposer4), TRUE, FALSE)) |> 
  select(Term, Bill_Vote_ID, BillTitle, leg_yea, URL, vote_pass, Proposer1 : Proposer9) |> 
  pivot_longer(Proposer1 : Proposer9) |> 
  filter(!is.na(value)) |> 
  select(Term, Bill_Vote_ID, BillTitle, vote_pass, leg_yea, Legislator = value, URL)
  
df_legislators_56 <- df_legislators |> 
  filter(term %in% c(5, 6)) |> 
  select(Legislator, election_type, Term = term)

df_bills_legis <- left_join(df_bills_legis, df_legislators_56)
```

Add bill names to rejected bills:

```{r}
df_bills_legis $ ShortBill <- NA
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_20_2013-11-07_2"] <- "Domestic Violence (2)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_21_2014-01-02_2"] <- "Rent Adjustment (2)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_22_2014-02-17_2"] <- "Animal Rights (2)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_23_2014-04-14_4"] <- "Trade Union Act (5)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_24_2014-07-01_6"] <- "Labour Protection (2)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_25_2014-08-13_11"] <- "Treaties on Workers' Rights (2)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_26_2015-06-15_9"] <- "Gender Equality in Employment (2)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_27_2015-06-15_10"] <- "Privacy Rights (2)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "2015/3_2015-01-13_5"] <- "* Organisation of the Legislature"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_28_2015-06-15_11"] <- "Labour Protection (3)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_29_2015-06-15_12"] <- "* Environmental Protection (2)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_30_2015-06-15_13"] <- "Treaties on Workers' Rights (3)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_31_2015-06-15_14"] <- "Trade Union Act (6)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "2017/13_2015-11-12_1"] <- "* Real Estate Leasing"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_32_2016-01-20_2"] <- "Trade Union Act (7)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_33_2016-01-21_3"] <- "* Penal Code (1)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_34_2016-11-10_9"] <- "Trade Union Act (8)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_35_2016-11-21_3"] <- "* Environmental Protection (3)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_36_2016-12-15_11"] <- "* Penal Code (2)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_37_2017-07-14_14"] <- "Treaties on Workers' Rights (4)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_38_2017-10-27_4"] <- "Trade Union Act (9)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_39_2019-04-23_3"] <- "Trade Union Act (10)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_40_2019-06-10_42"] <- "Treaties on Workers' Rights (5)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_41_2019-07-08_4"] <- "* Employment Policy"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_42_2020-03-16_23"] <- "Trade Union Act (11)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_43_2020-11-06_5"] <- "Trade Union Act (12)"
df_bills_legis $ ShortBill[df_bills_legis $ Bill_Vote_ID == "rejected_44_2021-04-08_18"] <- "Hearing of the Legislature"

df_bills_legis <- df_bills_legis |> 
  mutate(ShortBill = factor(ShortBill, levels = c(
    "* Real Estate Leasing", "* Organisation of the Legislature", 
    "Trade Union Act (7)", "Trade Union Act (8)", "* Penal Code (1)", 
    "* Penal Code (2)", "Animal Rights (2)",
    "Domestic Violence (2)", "Rent Adjustment (2)", "Trade Union Act (5)",
    "Trade Union Act (11)", "Trade Union Act (12)",
    "* Employment Policy", "* Environmental Protection (2)", "* Environmental Protection (3)",
    "Gender Equality in Employment (2)", "Hearing of the Legislature", 
    "Treaties on Workers' Rights (2)", "Treaties on Workers' Rights (3)", 
    "Treaties on Workers' Rights (4)", "Treaties on Workers' Rights (5)",
    "Privacy Rights (2)", "Trade Union Act (6)", "Trade Union Act (9)",
    "Trade Union Act (10)", "Labour Protection (2)", "Labour Protection (3)"
  ) ))

## Count the number of Yea's
df_bills_legis <- df_bills_legis |> 
  group_by(ShortBill) |> 
  mutate(n_sp = n())
```

## Plot

```{r}
#| label: figure-2
#| fig-cap: "Figure 2: Sponsors and supporters of legislator-proposed bills. This figure illustrates the sponsors of the 27 legislator-sponsored bills in the fifth and sixth terms. The bar height shows sponsors per bill; numbers above bars indicate “Yea” votes at the General Voting stage. Numbers in parentheses after a bill’s name denote how many similar bills have been submitted since the first term. An asterisk (*) indicates an amendment to existing law."
#| cap-location: margin
#| fig-width: 10
#| fig-height: 6

ggplot(data = df_bills_legis, 
       mapping = aes(x = ShortBill, fill = election_type)) + 
  geom_bar(stat = "count") + 
  geom_text(aes(y = n_sp, label = leg_yea, color = vote_pass), 
            vjust = 0, position = position_nudge(y = 0.5), show.legend = FALSE) + 
  scale_color_manual(values = c("darkgray", "black")) + 
  labs(x = NULL, y = "Number of sponsors", fill = "Legislator type") + 
  scale_y_continuous(breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9)) + 
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), 
        legend.position = c(0.9, 0.8),
        axis.text=element_text(size = 10),
        axis.title=element_text(size = 15),
        legend.title=element_text(size = 15),
        legend.text=element_text(size = 13)
        )
ggsave("Figures/Figure2_bill_proposers.pdf", width = 10, height = 6)
```

# Figure 3: Similarity Network for Term 5

Code trichotomous roll-call votes as ordinal:

```{r}
df_votes <- df_votes |> 
  mutate(vote_ord = ordered(vote, c("Yea", "Abstain", "Nay")))
with(df_votes, table(vote, vote_ord, useNA = "ifany"))
```

Recode the values so that absent = NA rather than 0.

```{r}
df_votes <- df_votes |> 
  mutate(vote_n = if_else(vote_n == 0, NA, vote_n))

with(df_votes, table(vote_n, vote_ord, useNA = "ifany"))
```

## Reshape

We reshape roll-call votes data from long to wide.

For Term 5:

```{r}
votes_wide_5 <- df_votes |> 
  filter(Term == 5) |>
  filter(Unanimous == FALSE) |> 
  select(Legislator, Bill_Vote_ID, vote_n) |>
  pivot_wider(names_from = "Bill_Vote_ID", 
              names_prefix = "bill_", 
              values_from = vote_n)  
```

This results in a 32 by 144 matrix.[^2]

[^2]: One bill was withdrawn after passing the General Voting stage.

For Term 6:

```{r}
votes_wide_6 <- df_votes |> 
  filter(Term == 6) |>
  filter(Unanimous == FALSE) |> 
  select(Legislator, Bill_Vote_ID, vote_n) |>
  pivot_wider(names_from = "Bill_Vote_ID", 
              names_prefix = "bill_", 
              values_from = vote_n)  
```

This results in a 33 by 208 matrix.

We add legislator information to these.

For Term 5:

```{r}
## Legislator information
df_legis_5 <- df_legislators |> 
  filter(term == 5) |>
  filter(Legislator != "賀一誠") |>
  select(Legislator, Legislator_E, LegislatorGender, election_type, election_type1,
         gender, party)

## Join
votes_wide_5 <- left_join(df_legis_5, votes_wide_5)
```

For Term 6:

```{r}
## Legislator information
df_legis_6 <- df_legislators |> 
  filter(term == 6) |>
  filter(Legislator != "賀一誠") |>
  select(Legislator, Legislator_E, LegislatorGender, election_type, election_type1,
         gender, party)

## Join
votes_wide_6 <- left_join(df_legis_6, votes_wide_6)
```

## Correlation Matrix

Calculate Spearman's rank correlation coefficients:

```{r}
cor_mat_pc <- votes_wide_5[-c(1:7)] |>
    t() |>
    cor(use = "pairwise.complete.obs", method = "spearman")

# Names
dimnames(cor_mat_pc) <- list(votes_wide_5 $ Legislator_E, votes_wide_5 $ Legislator_E)

# Fill the diagonal with 0s for the correlation of items with themselves
diag(cor_mat_pc) <- 0
```

Visualize the correlation matrix

```{r}
corrplot(cor_mat_pc, method = "circle", type = "full", tl.cex = 0.8, order = "AOE")
```

Network visualization

```{r}
cor_to_edge <- function(cormat, threshold = 0) {
  g <- graph_from_adjacency_matrix(cormat, mode = "upper", weighted = TRUE, diag = FALSE)
  e <- as_edgelist(g)
  df <- as.data.frame(cbind(e, E(g)$weight)) |>
    mutate(V3 = as.numeric(V3)) |>
    select(item1 = V1, item2 = V2, correlation = V3) |>
    filter(correlation > threshold)
}

graph_t5 <- cor_to_edge(cor_mat_pc, threshold = 0)

leglist_t5 <- graph_t5 |>
  select(item1, item2) |>
  unlist() |>
  unique()

parties <- votes_wide_5 |>
  select(Legislator = Legislator_E, LegislatorGender, election_type, election_type1,
         gender, party) |>
  filter(Legislator %in% leglist_t5 | Legislator %in% leglist_t5) |>
  mutate(party_o = recode(party, Peak = "1 Peak", Hometown = "2 Hometown", Gaming = "3 Gaming",
        `Pro-democracy` = "4 Pro-democracy", Other = "5 Other"))
```

## Plot

```{r}
#| label: figure-3
#| fig-cap: "Figure 3: Co-voting network for the fifth term (2013-2017). This figure presents a co-voting network for the fifth term. We assign different shapes to different types of political associations. Letters denote selection methods: “D” (directly elected), “I” (indirectly elected), and “A” (appointed). Edges indicate positive correlations, with thicker edges representing stronger similarities."
#| cap-location: margin
#| fig-width: 10
#| fig-height: 6

set.seed(1234567890)
graph_t5 |> 
  graph_from_data_frame(vertices = parties) |> 
  ggraph(layout = "stress") +
  geom_edge_link(aes(edge_alpha = correlation), color = "gray", show.legend = FALSE) +
  geom_node_point(aes(fill = party_o, color = party_o, shape = party_o), size = 5) +
  geom_node_text(aes(label = election_type1), color = "white", size = 4) +
  geom_node_text(aes(label = name), 
                 size = 3, check_overlap = F, vjust = 3, hjust = .4) +
  scale_shape_manual(values = c(15, 16, 17, 23, 25)) + 
  scale_color_viridis_d() + 
  scale_fill_viridis_d() + 
  coord_cartesian(clip = "off") + 
  scale_x_reverse() + 
  theme_void() + 
  labs(title = NULL, fill = "Association", 
       color = "Association", shape = "Association") + 
  theme(legend.position = "bottom",
        legend.key.size = unit(0.8, 'cm'), 
        legend.title = element_text(size = 17),
        legend.text = element_text(size = 15),
        plot.margin = unit(c(0, 30, 0, 30), "points")
        )

ggsave("Figures/Figure3_network_5.pdf", width = 10, height = 6)
```

## Distance

```{r}
leg_opp <- c("José Pereira Coutinho", 
             "Ng Kuok-cheong", 
             "Au Kam San", 
             "Leong Veng Chai", 
             "Kwan Tsui Hang", 
             "Lei Cheng I")
leg_est <- leglist_t5[leglist_t5 %in% leg_opp == FALSE]
```

### Within-camp correlations

Correlations between pro-establishment legislators

```{r}
graph_t5nf <- cor_to_edge(cor_mat_pc, threshold = -1)
graph_t5nf |> 
  filter(item1 %in% leg_est & item2 %in% leg_est) |> 
  pull(correlation) |> summary()
```

Correlations between anti-establishment legislators

```{r}
graph_t5nf |> 
  filter(item1 %in% leg_opp & item2 %in% leg_opp) |> 
  pull(correlation) |> summary()
```

### Between-camp correlations

Correlations between legislators from different groups

```{r}
graph_t5nf |> 
  filter(item1 %in% leg_opp & item2 %in% leg_est) |> 
  pull(correlation) |> summary()
```

# Figure 4: Similarity Network for Term 6

## Correlation Matrix

Calculate Spearman's rank correlation coefficients:

```{r}
cor_mat_pc <- votes_wide_6[-c(1:7)] |> 
  t() |> 
  cor(use = "pairwise.complete.obs", method = "spearman")

# Names
dimnames(cor_mat_pc) <- list(votes_wide_6 $ Legislator_E, votes_wide_6 $ Legislator_E)

# Fill the diagonal with 0s for the correlation of items with themselves
diag(cor_mat_pc) <- 0
```

Kou Hoi In and Wang Sai Man never overlap. Set their correlation as 0 instead of NA

```{r}
cor_mat_pc[rownames(cor_mat_pc) == "Kou Hoi In", colnames(cor_mat_pc) == "Wang Sai Man"] <- 0
cor_mat_pc[rownames(cor_mat_pc) == "Wang Sai Man", colnames(cor_mat_pc) == "Kou Hoi In"] <- 0
```

Visualize the correlation matrix

```{r}
corrplot(cor_mat_pc, method = "circle", type = "full", tl.cex = 0.8, order = "AOE")
```

Network visualization

```{r}
graph_t6 <- cor_to_edge(cor_mat_pc, threshold = 0)
leglist_t6 <- graph_t6 |> select(item1, item2) |> unlist() |> unique()
parties <- votes_wide_6 |> 
  select(Legislator = Legislator_E, LegislatorGender, election_type, election_type1, gender, party) |> 
  filter(Legislator %in% leglist_t6 | Legislator %in% leglist_t6) |> 
  mutate(party_o = recode(party, 
                          Peak = "1 Peak",
                          Hometown = "2 Hometown",
                          Gaming = "3 Gaming",
                          `Pro-democracy` = "4 Pro-democracy",
                          Other = "5 Other"
                          ))

```

## Plot

```{r}
#| label: figure-4
#| fig-cap: "Figure 4: Co-voting network for the sixth term (2017-2021). This figure presents a co-voting network for the sixth term. Please refer to Figure 3 for details."
#| cap-location: margin
#| fig-width: 10
#| fig-height: 6

set.seed(1234567890)
graph_t6 |> 
  graph_from_data_frame(vertices = parties) |> 
  ggraph(layout = "stress") +
  geom_edge_link(aes(edge_alpha = correlation), color = "gray", show.legend = FALSE) +
  geom_node_point(aes(fill = party_o, color = party_o, shape = party_o), size = 5) +
  geom_node_text(aes(label = election_type1), color = "white", size = 4) +
  geom_node_text(aes(label = name), 
                 size = 3, check_overlap = F, vjust = 3, hjust = .4) +
  scale_shape_manual(values = c(15, 16, 17, 23, 25)) + 
  scale_color_viridis_d() + 
  scale_fill_viridis_d() + 
  coord_cartesian(clip = "off") + 
  theme_void() + 
  labs(title = NULL, fill = "Association", 
       color = "Association", shape = "Association") + 
  theme(legend.position = "bottom",
        legend.key.size = unit(0.8, 'cm'), 
        legend.title = element_text(size = 17),
        legend.text = element_text(size = 15),
        plot.margin = unit(c(0, 30, 0, 30), "points")
        )
ggsave("Figures/Figure4_network_6.pdf", width = 10, height = 6)
```

## Distance

```{r}
leg_opp <- c("José Pereira Coutinho", 
             "Sou Ka Hou",
             "Ng Kuok-cheong", 
             "Au Kam San")
leg_est <- leglist_t6[leglist_t6 %in% leg_opp == FALSE]
```

### Within-camp correlations

Correlations between pro-establishment legislators

```{r}
graph_t6nf <- cor_to_edge(cor_mat_pc, threshold = -1)
graph_t6nf |> 
  filter(item1 %in% leg_est & item2 %in% leg_est) |> 
  pull(correlation) |> summary()
```

Correlations between anti-establishment legislators

```{r}
graph_t6nf |> 
  filter(item1 %in% leg_opp & item2 %in% leg_opp) |> 
  pull(correlation) |> summary()
```

### Between-camp correlations

Correlations between legislators from different groups

```{r}
graph_t6nf |> 
  filter(item1 %in% leg_opp & item2 %in% leg_est) |> 
  pull(correlation) |> summary()
```

# Figure 5: Ideal Points for Term 5

Roll-call votes data for Term 5:

```{r}
df_votes5 <- df_votes |> 
  filter(Unanimous == FALSE) |> 
  filter(Term == 5) |> 
  filter(!is.na(vote)) |> 
  select(Legislator, vote, vote_ord, Bill_Vote_ID)
```

## Bayesian IRT Model

Set up a two parameter ordinal IRT model:

```{r}
## Formula
formula_va_ord_2pl <- bf(
  vote_ord ~ alpha * eta, 
  eta ~ 1 + (1 |i| Bill_Vote_ID) + (1 | Legislator),
  alpha ~ 1 + (1 |i| Bill_Vote_ID),
  nl = TRUE
)

## Prior: sd for theta = 1
prior_va_ord_2pl <- 
  set_prior("constant(1)", class = "sd", group = "Legislator", nlpar = "eta") + 
  set_prior("normal(0, 1)", class = "sd", group = "Bill_Vote_ID", nlpar = "eta") + 
  set_prior("normal(0, 1)", class = "sd", group = "Bill_Vote_ID", nlpar = "alpha")
```

`(1 |i| Bill_Vote_ID)` allows for bill-level parameters (difficulty and discrimination) to be correlated. For example, more difficult parameters may have higher discrimination.

The following code chunk estimates the model for the fifth term. On my machine (a 2024 MacBook Pro with an Apple M4 processor), this process took about 30 minutes. I have set the code chunk option to `eval = FALSE` so that the code will not run by default. Please set it to `eval = TRUE` if you wish to execute the code.

```{r eval = FALSE}
fit_term5 <- brm(
  formula = formula_va_ord_2pl,
  data = df_votes5,
  family = brmsfamily("cumulative", "logit"),
  prior = prior_va_ord_2pl,
  chains = 1,
  iter = 10000,
  cores = 10,
  seed = 12345,
  control = list(adapt_delta = 0.95, max_treedepth = 15)
)
write_rds(fit_term5, file = "Data/fit_term5_10k.RDS")
```

Load the estimation results:

```{r}
fit_term5 <- read_rds("Data/fit_term5_10k.RDS")
summary(fit_term5)
```

## Plot Legislators Ideal Points

Legislator's election type, gender, and party

```{r}
coef_eta <- coef(fit_term5) $ Legislator[,,1] |> 
  as.data.frame() |> 
  mutate(Legislator = rownames(coef(fit_term5) $ Legislator[,,1])) |> 
  mutate(Mid = Estimate, Upper = Q97.5, Lower = Q2.5)

## Merge legislator names
coef_eta <- left_join(df_legis_5, coef_eta)

## Adjust x-axis position on labels
coef_eta <- coef_eta |> 
  mutate(xpos = if_else(Lower > 0, Upper - 2, Lower + 2.5))
## Adjust name location for those with long names
coef_eta $ xpos[coef_eta $ Legislator == "高天賜"] <- coef_eta $ xpos[coef_eta $ Legislator == "高天賜"] - 0.5
```

```{r}
#| label: figure-5
#| fig-cap: "Figure 5: Estimated ideal points of legislators for the fifth legislative term, 2013-2017. This figure presents ideal points estimated from roll-call votes for the fifth term. We assign different shapes to different types of political associations. Letters denote selection methods: “D” (directly elected), “I” (indirectly elected), and “A” (appointed). Lower values indicate pro-establishment positions; higher values indicate anti-establishment positions. Horizontal lines show 95% Bayesian credible intervals. Women legislators are denoted with an asterisk (*) at the end of their name."
#| cap-location: margin
#| fig-width: 10
#| fig-height: 6

ggplot(data = coef_eta,
        mapping = aes(y = reorder(LegislatorGender, Mid),
                      x = Mid, color = party)) +
   geom_segment(aes(x = Lower, xend = Upper,
                    yend = reorder(LegislatorGender, Mid)),
                linewidth = 1.5, show.legend = FALSE) +
   geom_text(aes(label = LegislatorGender, x = xpos), size = 3.5, color = "gray40") +
   geom_point(aes(shape = party, fill = party), size = 4) +
   geom_text(aes(label = election_type1), color = "white", size = 3, show.legend = FALSE) +
   scale_shape_manual(values = c(15, 16, 17, 23, 25)) +
   scale_color_viridis_d() +
   scale_fill_viridis_d() +
   labs(x = "Ideal points",
        y = NULL,
        # title = "Estimated ideal points of Macao SAR legislators",
        # subtitle = "Term 5 (2013 - 2017)",
        color = "Association",
        shape = "Association",
        fill = "Association"
        ) +
   theme_minimal() +
   theme(legend.position = "bottom",
         axis.text.y = element_blank(),
         legend.key.size = unit(0.8, 'cm'),
         strip.text = element_text(size = 15),
         axis.text.x = element_text(size = 10),
         axis.title = element_text(size = 17),
         legend.title = element_text(size = 17),
         legend.text = element_text(size = 15)
         )

ggsave("Figures/Figure5_idealpoints_5.pdf", width = 10, height = 6)
```

# Figure 6: Ideal Points for Term 6

Roll-call votes data for Term 6:

```{r}
df_votes6 <- df_votes |> 
  filter(Unanimous == FALSE) |> 
  filter(Term == 6) |> 
  filter(!is.na(vote)) |> 
  select(Legislator, vote, vote_ord, Bill_Vote_ID)
```

## Bayesian IRT Model

The following code chunk estimates the model for the fifth term. On my machine (a 2024 MacBook Pro with an Apple M4 processor), this process took about 240 minutes. I have set the code chunk option to `eval = FALSE` so that the code will not run by default. Please set it to `eval = TRUE` if you wish to execute the code.

I increased the number of iterations to 20,000 for the sixth term to improve model fit, as indicated by the fit statistics.

```{r eval = FALSE}
fit_term6 <- brm(
  formula = formula_va_ord_2pl,
  data = df_votes6,
  family = brmsfamily("cumulative", "logit"),
  prior = prior_va_ord_2pl,
  chains = 1,
  iter = 20000,
  cores = 10,
  seed = 12345,
  control = list(adapt_delta = 0.95, max_treedepth = 15)
)

write_rds(fit_term6, file = "Data/fit_term6_20k.RDS")
```

```{r}
fit_term6 <- read_rds("Data/fit_term6_20k.RDS")
summary(fit_term6)
```

## Plot Legislators Ideal Points

Legislator's election type, gender, and party

```{r}
coef_eta <- coef(fit_term6) $ Legislator[,,1] |> as.data.frame() |> 
  mutate(Legislator = rownames(coef(fit_term6) $ Legislator[,,1])) |> 
  mutate(Mid = Estimate, Upper = Q97.5, Lower = Q2.5)

## Merge legislator names
coef_eta <- left_join(df_legis_6, coef_eta)

coef_eta <- coef_eta |> 
  mutate(xpos = if_else(Mid > 0, Lower - 1.5, Upper + 1))
```

```{r}
#| label: figure-6
#| fig-cap: "Figure 6: Estimated ideal points of legislators for the sixth legislative term, 2017-2021. This figure presents legislators’ ideal points for the sixth term, with 95% Bayesian credible intervals. See Figure 5 for details."
#| cap-location: margin
#| fig-width: 10
#| fig-height: 6

ggplot(data = coef_eta, 
       mapping = aes(y = reorder(LegislatorGender, Mid), 
                     x = Mid, color = party)) + 
  geom_segment(aes(x = Lower, xend = Upper, 
                   yend = reorder(LegislatorGender, Mid)), 
               linewidth = 1.5, show.legend = FALSE) + 
  geom_text(aes(label = LegislatorGender, x = xpos), size = 3.5, color = "gray40") + 
  geom_point(aes(shape = party, fill = party), size = 4) + 
  geom_text(aes(label = election_type1), color = "white", size = 3, show.legend = FALSE) + 
  scale_shape_manual(values = c(15, 16, 17, 23, 25)) + 
  scale_color_viridis_d() + 
  scale_fill_viridis_d() + 
  labs(x = "Ideal points", 
       y = NULL, 
       color = "Association",
       shape = "Association",
       fill = "Association"
       ) + 
  theme_minimal() + 
  theme(legend.position = "bottom",
        axis.text.y = element_blank(),
        legend.key.size = unit(0.8, 'cm'), 
        strip.text = element_text(size = 15),
        axis.text.x = element_text(size = 10),
        axis.title = element_text(size = 17),
        legend.title = element_text(size = 17),
        legend.text = element_text(size = 15)
        )

ggsave("Figures/Figure6_idealpoints_6.pdf", width = 10, height = 6)
```

# Figure 7: Deliberation Duration

## Data Cleaning

Our analysis will include those bills that:

-   were proposed by the executive government (= we exclude legislator-proposed bills);
-   passed the General Voting stage; and
-   were assigned to one of the three committees.

Collect names of rejected bills and withdrawn bills:

```{r}
## Rejected
rej_bills <- c()
for(i in 1:44){
  rej_bills <- c(rej_bills, paste("rejected_", i, sep = ""))
}

## Withdrawn
withdrawn_bills <- df_bills $ BillID[grep("政府撤回", df_bills $ BillID)]
```

```{r}
df_bills_gd <- df_bills |> 
  filter(BillID %in% rej_bills == FALSE) |> 
  filter(BillID %in% withdrawn_bills == FALSE)

nrow(df_bills_gd)
```

Split General and Detailed:

```{r}
df_bills_gen <- df_bills_gd |> 
  filter(vote_type == "General") |> 
  filter(proposer == "Government")

df_bills_det <- df_bills_gd |> 
  filter(vote_type == "Detailed")

with(df_bills_gen, table(vote_unanimous, Term, useNA = "ifany"))
with(df_bills_det, table(vote_unanimous, Term, useNA = "ifany"))
```

## Covariates

General voting with oppositions:

```{r}
dfg56 <- df_bills_gen |> 
  filter(vote_unanimous == FALSE & gov_bill == "Yes" & com_assign != "Not Assigned") |> 
  select(Bill_Vote_ID, Term)
```

Nay & abstain votes:

```{r}
dfv <- df_votes |> 
  filter(Unanimous == FALSE & vote_type == "一般性討論") |> 
  filter(Bill_Vote_ID %in% dfg56 $ Bill_Vote_ID) |> 
  filter(vote_ord != "Yea" & !is.na(vote_ord)) |> 
  select(Legislator, Term, Bill_Vote_ID, vote_ord)
```

Extract bills that had opposition votes (nay or abstain) from pro-establishment legislators:

```{r}
dfv_opp <- dfv |> 
  filter(Legislator != "吳國昌") |> 
  filter(Legislator != "區錦新") |> 
  filter(Legislator != "高天賜") |> 
  filter(Legislator != "梁榮仔") |> 
  filter(Legislator != "蘇嘉豪") |> 
  filter(Legislator != "李靜儀" | Term == 6) |> 
  filter(Legislator != "關翠杏") |> 
  pull(Bill_Vote_ID) |> unique()
```

Choose one detailed vote per bill:

```{r}
df_bills_det1 <- df_bills_det |> 
  group_by(BillID) |> 
  summarise(vote_date_detail = max(vote_date))
```

Merge General- and Detailed-stage votes:

```{r}
df_bills_dur <- tidylog::left_join(df_bills_gen, df_bills_det1)
```

The two unmatched cases are 2015/3 and 2017/13, both of which are legislator proposed.

### Opposition by estalishment (opp_est)

```{r}
df_bills_dur <- df_bills_dur |> 
  mutate(opp_est = if_else(Bill_Vote_ID %in% dfv_opp, "Yes", "No"))

with(df_bills_dur, table(opp_est, useNA = "ifany"))
```

Check the numbers

```{r}
with(df_bills_dur, table(vote_unanimous, Term, useNA = "ifany"))
```

### Any Opposition

```{r}
df_bills_dur <- df_bills_dur |> 
  mutate(opp_any = if_else(vote_unanimous == FALSE, "Yes", "No"), 
         opp_opp = if_else(opp_any == "Yes" & opp_est == "No", "Yes", "No"))

with(df_bills_dur, table(opp_est, opp_any, useNA = "ifany"))
```

### Amendment Bills

```{r}
bill_am <- df_bills_dur $ BillTitle[grep("修改", df_bills_dur $ BillTitle)]
head(bill_am)
```

```{r}
df_bills_dur <- df_bills_dur |> 
  mutate(am_bill = if_else(BillTitle %in% bill_am, "Yes", "No"))
table(df_bills_dur $ am_bill)
```

## Bills Not Assigned to Any Committee

Show those bills that are not assigned to a committee:

```{r}
df_bills_dur |> 
  filter(com_assign %in% c("1", "2", "3") == FALSE) |> 
  kable()
```

Among the 14 (out of 121 gov-proposed bills),

-   4 involve amendment to 禁止不法生産 bills.
-   3 involve amendment to annual budget.
-   3 involve amendment to 規章 of 消費稅, 機動車輛稅, 市區房屋稅
-   1 involves 抗日戰爭 anniversary.
-   1 on 取得非首個居住用途不動產的印花稅.
-   1 involves 附件的傳染病表.
-   1 involves 批准澳門特別行政區政府承擔

```{r}
table(df_bills_dur $ com_assign, useNA = "ifany")
```

We drop those bills not assigned to committee.

```{r}
df_bills_dur <- df_bills_dur |> 
  filter(com_assign %in% c("1", "2", "3"))

dim(df_bills_dur)
```

Clean-up the data:

```{r}
df_tc <- df_bills_dur |> 
  mutate(duration = vote_date_detail - vote_date + 1, 
         duration_num = as.numeric(duration)) |> 
  select(-c(Proposer1:URL))
```

## Descriptive Statistics

```{r}
ggplot(data = df_tc, 
       mapping = aes(x = as.factor(Term), 
                     y = duration_num, 
                     fill = vote_unanimous)) + 
   geom_boxplot() + 
  labs(x = "Term", y = "Duration (days)")
```

```{r}
with(df_tc, by(duration_num, Term, summary))
```

```{r}
summary(df_tc $ duration_num)
```

```{r}
df_tc |> filter(duration_num < 180) |> nrow()
```

63 out of 121 (52%) of bills finished in less than six months.

```{r}
df_tc |> filter(duration_num > 365) |> nrow()
27/121
```

27 out of 121 (22%) of bills took more than a year.

## Plot

```{r}
#| label: figure-7
#| fig-cap: "Figure 7: Committee deliberation duration for government-sponsored bills. This figure shows the distribution of committee deliberation duration in days, measured as the time between General Voting and Detailed Voting stages."
#| cap-location: margin
#| fig-width: 10
#| fig-height: 6

ggplot(data = df_tc, 
       mapping = aes(x = duration_num)) + 
  geom_histogram(color = "black", 
                 fill = viridis::viridis(1), 
                 boundary = 0.5) + 
  labs(x = "Deliberation duration (days)", 
       y = "Number of bills") + 
  theme_classic() + 
  theme(axis.text=element_text(size = 15),
        axis.title=element_text(size = 17))
ggsave("Figures/Figure7_duration.pdf", width = 10, height = 6)
```

# Table 2: Survival Analysis

We perform survival analysis using Cox duration models. As our analysis incorporates time-varying covariates, we convert the duration data into the time-varying format.

## Convert to Time-Varying Data

```{r}
df_tv <- df_tc |> 
  rowwise() |> 
  mutate(dates = list(seq.Date(vote_date, vote_date_detail, by = "day"))) |> 
  unnest(dates) |> 
  ungroup() |> 
  mutate(day = as.numeric(dates - vote_date + 1))
```

### Covariates

We define the following covariates:

1.  Number of days left in session
2.  Number of bills in total at a given time
3.  Number of bills in the same committee at a given time

#### Days left

Session is from October 16 to October 15 of the following year, so

-   Term 5 is from 16 October 2013 to 15 October 2017.
-   Term 6 is from 16 October 2017 to 15 October 2021.

```{r}
table(df_tv $ Term)
```

```{r}
df_tv <- df_tv |> 
  mutate(days_left = if_else(Term == 5, 
                             as.Date("2017-10-15") - dates,
                             as.Date("2021-10-15") - dates
                             ), 
         days_left_num = as.numeric(days_left), 
         days_left_log = log(days_left_num), 
         NumPages_log = log(NumPages))
```

#### Number of bills in total

Count the number of **other** bills that are active on a given day.

```{r}
df_tv <- df_tv |> 
  group_by(dates) |> 
  mutate(total_bill_day = n() - 1) |> 
  ungroup()

ggplot(data = df_tv, aes(x = total_bill_day)) + geom_histogram()
```

#### Number of bills in committee

Count the number of **other** bills in the same committee that are active on a given day.

```{r}
df_tv <- df_tv |> 
  group_by(dates, com_assign) |> 
  mutate(com_bill_day = n() - 1) |> 
  ungroup()

table(df_tv $ com_bill_day)
```

## Survival Analysis

Define some helper functions

```{r, echo = FALSE}
## Functions to calculate clustered standard errors from the fitted object
### For hazard rate
get_clse_hr <- function(fit){
  sqrt(diag(fit $ var))
  }

### For hazard ratio (exponentiated coefficients)
get_clse_or <- function(fit){
  or <- exp(coef(fit))
  var.diag <- diag(fit $ var)
  return(sqrt(or^2 * var.diag))
}

# ## Functions to extract relevant results
# extract_hr <- function(fit){
#   df <- summary(fit) $ conf.int
#   pe <- df[rownames(df) == "crisis_mon", "exp(coef)"]
#   low <- df[rownames(df) == "crisis_mon", "lower .95"]
#   up <- df[rownames(df) == "crisis_mon", "upper .95"]
#   data.frame(pe = pe, low = low, up = up)
# }

## Function to add AIC to the model object (so that stargazer recognizes them)
add_AIC <- function(fit){
  fit $ AIC <- AIC(fit)
  return(fit)
}

zphTable <- function(cox.fit, digits) {
  # Performs cox.zph function from the survival package using the four
  #   available time transformations and outputs one table with side-by-side
  #   presentation of results
  # Args:
  #   cox.fit: fit of a Cox model using the coxph function
  #   digits:  desired number of digits for presentation of results
  tab <- round(cbind(cox.zph(cox.fit, transform = 'km')$table,
                     cox.zph(cox.fit, transform = 'rank')$table), digits = digits)
  cols <- colnames(cox.zph(cox.fit)$table)
  colnames(tab) <- c(paste('km.', cols, sep = ''),
                     paste('rank.', cols, sep = ''))    
  tab
}
```

Define event indicator: 1 = event occurred, 0 = censored

```{r}
## Time-constant data
df_tc $ event <- 1
## Time-varying data
df_tv <- df_tv |> 
  mutate(event = if_else(dates == vote_date_detail, 1, 0))
# table(df_tc $ event)
# table(df_tv $ event)
```

Generate `day0`, which marks the beginning (t - 1).

```{r}
df_tv <- df_tv |> 
  mutate(day0 = day - 1)
```

## Cox Models

### Compare Time-Constant and Time-Varying Models

We first compare the time-constant and time-varying results (they should be identical if the data conversion was correctly done).

```{r}
## Time-constant
srv_tc <- with(df_tc, Surv(duration, event = event))
## Time-varying
srv_tv <- with(df_tv, Surv(time = day0, time2 = day, event = event))

fit_tc1 <- coxph(srv_tc ~ am_bill + opp_est + as.factor(Term), 
                 data = df_tc, cluster = BillID)
fit_tc2 <- coxph(srv_tc ~ am_bill + opp_opp + as.factor(Term), 
                 data = df_tc, cluster = BillID)
fit_tv1 <- coxph(srv_tv ~ am_bill + opp_est + as.factor(Term), 
                 data = df_tv, cluster = BillID)
fit_tv2 <- coxph(srv_tv ~ am_bill + opp_opp + as.factor(Term), 
                 data = df_tv, cluster = BillID)

model_list <- list(fit_tc1, fit_tc2, fit_tv1, fit_tv2)
## Add AIC
model_list <- lapply(model_list, add_AIC)
```

```{r}
#| warning: false
#| results: asis

stargazer(model_list,
          se = lapply(model_list, get_clse_hr),
          omit.stat = c("rsq", "max.rsq", "wald", "lr", "logrank"),
          type = "html")
```

Time-constant results (models 1 and 2) are indeed identical with time-varying results (models 3 and 4).

### Opposition by pro-establishment

```{r}
srv_tv5 <- with(df_tv |> filter(Term == 5), 
                Surv(time = day0, time2 = day, event = event))
srv_tv6 <- with(df_tv |> filter(Term == 6), 
                Surv(time = day0, time2 = day, event = event))
df_tv <- df_tv |> 
  mutate(term6 = if_else(Term == 6, 1, 0))
```

```{r}
fit_tv1 <- coxph(srv_tv ~ opp_est + term6 + 
                   NumPages_log + days_left_log + com_bill_day, 
                 data = df_tv, 
                 cluster = BillID)
fit_tv1_5 <- coxph(srv_tv5 ~ opp_est + 
                     NumPages_log + days_left_log + com_bill_day, 
                 data = df_tv |> filter(Term == 5), 
                 cluster = BillID)
fit_tv1_6 <- coxph(srv_tv6 ~ opp_est + 
                   NumPages_log + days_left_log + com_bill_day, 
                 data = df_tv |> filter(Term == 6), 
                 cluster = BillID)

model_list <- list(fit_tv1, fit_tv1_5, fit_tv1_6)
## Add AIC
model_list <- lapply(model_list, add_AIC)
```

```{r, echo = FALSE}
#| warning: false
#| results: asis
stargazer(model_list,
          se = lapply(model_list, get_clse_hr),
          omit.stat = c("rsq", "max.rsq", "wald", "lr", "logrank"),
          type = "html")
```

#### Testing non-proportionality

We follow Park and Hendry’s (2015) approach to test the non-proportionality using different transformations of time.

```{r}
zphTable(fit_tv1, digits = 3)
```

```{r}
zphTable(fit_tv1_5, digits = 3)
```

```{r}
zphTable(fit_tv1_6, digits = 3)
```

Include time interactions

```{r}
fit_ti1 <- coxph(srv_tv ~ opp_est + 
                   NumPages_log + com_bill_day + days_left_log + 
                   term6 + tt(term6), 
                 tt = function(x, t, ...) x * log(t),
                 data = df_tv, 
                 cluster = BillID)

fit_ti1_5 <- coxph(srv_tv5 ~ opp_est + 
                     NumPages_log + com_bill_day + 
                     days_left_log + tt(days_left_log), 
                 tt = function(x, t, ...) x * log(t),
                   data = df_tv |> filter(Term == 5), 
                   cluster = BillID)
fit_ti1_6 <- coxph(srv_tv6 ~ opp_est + 
                     NumPages_log + com_bill_day + tt(com_bill_day) + days_left_log, 
                 tt = function(x, t, ...) x * log(t),
                   data = df_tv |> filter(Term == 6), 
                   cluster = BillID)

model_list <- list(fit_ti1, fit_ti1_5, fit_ti1_6)
## Add AIC
model_list <- lapply(model_list, add_AIC)
```

```{r, echo = FALSE}
#| warning: false
#| results: asis
stargazer(model_list,
          se = lapply(model_list, get_clse_hr),
          omit.stat = c("rsq", "max.rsq", "wald", "lr", "logrank"),
          type = "html")
```

#### Opposition by opposition

```{r}
fit_tv2 <- coxph(srv_tv ~ opp_opp + term6 + 
                   NumPages_log + days_left_log + com_bill_day, 
                 data = df_tv, 
                 cluster = BillID)
fit_tv2_5 <- coxph(srv_tv5 ~ opp_opp + 
                     NumPages_log + days_left_log + com_bill_day, 
                 data = df_tv |> filter(Term == 5), 
                 cluster = BillID)
fit_tv2_6 <- coxph(srv_tv6 ~ opp_opp + 
                   NumPages_log + days_left_log + com_bill_day, 
                 data = df_tv |> filter(Term == 6), 
                 cluster = BillID)

model_list <- list(fit_tv2, fit_tv2_5, fit_tv2_6)
## Add AIC
model_list <- lapply(model_list, add_AIC)
```

```{r, echo = FALSE}
#| warning: false
#| results: asis
stargazer(model_list,
          se = lapply(model_list, get_clse_hr),
          omit.stat = c("rsq", "max.rsq", "wald", "lr", "logrank"),
          type = "html")
```

##### Testing non-proportionality

We follow Park and Hendry's (2015) approach to test the non-proportionality using different transformations of time.

```{r}
zphTable(fit_tv2, digits = 3)
```

Model 1 (term 5)

```{r}
zphTable(fit_tv2_5, digits = 3)
```

Model 1 (term 6)

```{r}
zphTable(fit_tv2_6, digits = 3)
```

##### Time interactions

```{r}
fit_ti2 <- coxph(srv_tv ~ opp_opp + 
                   NumPages_log + com_bill_day + 
                   days_left_log + term6 + tt(term6), 
                 tt = function(x, t, ...) x * log(t),
                 data = df_tv, 
                 cluster = BillID)

fit_ti2_5 <- coxph(srv_tv5 ~ opp_opp + 
                     NumPages_log + com_bill_day + 
                     days_left_log + tt(days_left_log), 
                 tt = function(x, t, ...) x * log(t),
                   data = df_tv |> filter(Term == 5), 
                   cluster = BillID)
fit_ti2_6 <- coxph(srv_tv6 ~ opp_opp + 
                     NumPages_log + 
                     com_bill_day + tt(com_bill_day) + 
                     days_left_log, 
                 tt = function(x, t, ...) x * log(t),
                   data = df_tv |> filter(Term == 6), 
                   cluster = BillID)

model_list <- list(fit_ti2, fit_ti2_5, fit_ti2_6)
## Add AIC
model_list <- lapply(model_list, add_AIC)
```

```{r, echo = FALSE}
#| warning: false
#| results: asis
stargazer(model_list,
          se = lapply(model_list, get_clse_hr),
          omit.stat = c("rsq", "max.rsq", "wald", "lr", "logrank"),
          type = "html")
```

#### Opposition by any

```{r}
fit_tv3 <- coxph(srv_tv ~ opp_any + term6 + 
                   NumPages_log + days_left_log + com_bill_day, 
                 data = df_tv, 
                 cluster = BillID)
fit_tv3_5 <- coxph(srv_tv5 ~ opp_any + 
                     NumPages_log + days_left_log + com_bill_day, 
                 data = df_tv |> filter(Term == 5), 
                 cluster = BillID)
fit_tv3_6 <- coxph(srv_tv6 ~ opp_any + 
                   NumPages_log + days_left_log + com_bill_day, 
                 data = df_tv |> filter(Term == 6), 
                 cluster = BillID)

model_list <- list(fit_tv3, fit_tv3_5, fit_tv3_6)
## Add AIC
model_list <- lapply(model_list, add_AIC)
```

```{r, echo = FALSE}
#| warning: false
#| results: asis
stargazer(model_list,
          se = lapply(model_list, get_clse_hr),
          omit.stat = c("rsq", "max.rsq", "wald", "lr", "logrank"),
          type = "html")
```

##### Testing non-proportionality

We follow Park and Hendry's (2015) approach to test the non-proportionality using different transformations of time.

```{r}
zphTable(fit_tv3, digits = 3)
```

Model 1 (term 5)

```{r}
zphTable(fit_tv3_5, digits = 3)
```

Model 1 (term 6)

```{r}
zphTable(fit_tv3_6, digits = 3)
```

##### Time interactions

```{r}
fit_ti3 <- coxph(srv_tv ~ opp_any + 
                   NumPages_log + com_bill_day + days_left_log + 
                   term6 + tt(term6), 
                 tt = function(x, t, ...) x * log(t),
                 data = df_tv, 
                 cluster = BillID)

fit_ti3_5 <- coxph(srv_tv5 ~ opp_any + 
                     NumPages_log + 
                     com_bill_day + tt(com_bill_day) + 
                     days_left_log, 
                 tt = function(x, t, ...) x * log(t),
                   data = df_tv |> filter(Term == 5), 
                   cluster = BillID)
fit_ti3_6 <- coxph(srv_tv6 ~ opp_any + 
                     NumPages_log + com_bill_day + 
                     days_left_log + tt(days_left_log), 
                 tt = function(x, t, ...) x * log(t),
                   data = df_tv |> filter(Term == 6), 
                   cluster = BillID)

model_list <- list(fit_ti3, fit_ti3_5, fit_ti3_6)
## Add AIC
model_list <- lapply(model_list, add_AIC)
```

```{r, echo = FALSE}
#| warning: false
#| results: asis
stargazer(model_list,
          se = lapply(model_list, get_clse_hr),
          omit.stat = c("rsq", "max.rsq", "wald", "lr", "logrank"),
          type = "html")
```

## Results Summary

```{r}
fit_tv4 <- coxph(srv_tv ~ opp_est + opp_opp + term6 + 
                   NumPages_log + days_left_log + com_bill_day, 
                 data = df_tv, 
                 cluster = BillID)

```

```{r, echo = FALSE, eval = FALSE, warning = FALSE}
model_list <- list(fit_tv1, fit_tv2, fit_tv4)
## Add AIC
model_list <- lapply(model_list, add_AIC)
stargazer(model_list,
          se = lapply(model_list, get_clse_hr),
          omit.stat = c("rsq", "max.rsq", "wald", "lr", "logrank"),
          covariate.labels = c("Opposition by pro-establishment", 
                               "Opposition by opposition",
                               "Term 6",
                               "Bill Length (logged)", 
                               "Days left in session (logged)", 
                               "Number of bills"),
          dep.var.caption = c("Hazard Rate (> 0 means shorter duration)"),
          dep.var.labels.include = FALSE,
          title = "Deliberation",
          type = "text")
```

```{r}
exp(-1.143)
```

Opposition by a pro-establishment legislator reduces the hazard by a factor of 0.32, or 68%.

```{r, warning = FALSE, results = 'asis'}
model_list <- list(fit_tv1, fit_tv2, fit_tv4)
## Add AIC
model_list <- lapply(model_list, add_AIC)
stargazer(model_list,
          se = lapply(model_list, get_clse_hr),
          omit.stat = c("rsq", "max.rsq", "wald", "lr", "logrank"),
          covariate.labels = c("Opposition by pro-establishment", 
                               "Opposition by opposition",
                               "Term 6",
                               "Bill Length (logged)", 
                               "Days left in session (logged)", 
                               "Number of bills"),
          dep.var.caption = c("Hazard Rate (> 0 means shorter duration)"),
          dep.var.labels.include = FALSE,
          title = "Table 2: Committee Deliberation Duration",
          type = "html")
```

```{r}
fit_tv4_5 <- coxph(srv_tv5 ~ opp_est + opp_opp + 
                   NumPages_log + days_left_log + com_bill_day, 
                 data = df_tv |> filter(Term == 5), 
                 cluster = BillID)
fit_tv4_6 <- coxph(srv_tv6 ~ opp_est + opp_opp + 
                   NumPages_log + days_left_log + com_bill_day, 
                 data = df_tv |> filter(Term == 6), 
                 cluster = BillID)

```

```{r, echo = FALSE, eval = FALSE, warning = FALSE}
model_list <- list(fit_tv1_5, fit_tv2_5, fit_tv4_5, fit_tv1_6, fit_tv2_6, fit_tv4_6)
## Add AIC
model_list <- lapply(model_list, add_AIC)
stargazer(model_list,
          se = lapply(model_list, get_clse_hr),
          omit.stat = c("rsq", "max.rsq", "wald", "lr", "logrank"),
          covariate.labels = c("Opposition by pro-establishment", 
                               "Opposition by opposition",
                               "Bill Length (logged)", 
                               "Days left in session (logged)", 
                               "Number of bills"),
          dep.var.caption = c("Hazard Rate (> 0 means shorter duration)"),
          dep.var.labels.include = FALSE,
          title = "Robustness check: separate analyses for terms",
          type = "text")
```

```{r}
#| results: asis
stargazer(model_list,
          se = lapply(model_list, get_clse_hr),
          omit.stat = c("rsq", "max.rsq", "wald", "lr", "logrank"),
          covariate.labels = c("Opposition by pro-establishment", 
                               "Opposition by opposition",
                               "Bill Length (logged)", 
                               "Days left in session (logged)", 
                               "Number of bills"),
          dep.var.caption = c("Hazard Rate (> 0 means shorter duration)"),
          column.labels = c("Term 5", "Term 5", "Term 5", "Term 6", "Term 6", "Term 6"),
          dep.var.labels.include = FALSE,
          title = "Robustness check: separate analyses for terms",
          type = "html")
```

<details>

<summary>Session Info</summary>

```{r sessionInfo}
Sys.time()
sessionInfo()
```

</details>
